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\ 1 INTRODUCTION 



ABSTRACT 

We use three-dimensional SPH simulations to investigate the collapse of low-mass 
prestellar cores and the formation and early evolution of protostellar discs. The initial 
conditions are slightly supercritical Bonnor-Ebert spheres in rigid rotation. The core 
mass and initial radius are held fixed at M Q = 6.1 Mq and R Q = 17, 000 AU, and the 
only parameter that we vary is the initial angular speed il Q . Protostellar discs forming 
from cores with fl Q < 1.35 x 10~ 13 s _1 have radii between 100 and 300 AU and are 
quite centrally concentrated; due to heating by gas infall onto the disc and accretion 
onto the central object, they are also quite warm, T > 100 K, and therefore stable 
against gravitational fragmentation. In contrast, more rapidly rotating cores form discs 
which are less concentrated and cooler, and have radii between 400 and 1000 AU; as 
a consequence they are prone to gravitational fragmentation and the formation of 
multiple systems. We derive a criterion that predicts whether a rigidly rotating core 
having given M Q , R Q and f2 Q will produce a protostellar disc which fragments whilst 
material is still infalling from the core envelope. We then apply this criterion to core 
samples for which M Q , R Q and o have been estimated observationally. We conclude 
that the observed cores are stable against fragmentation at this stage, due to their 
low angular speeds and the heat delivered at the accretion shock where the infalling 
material hits the disc. 
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Stars form from the gravitational collapse of prestellar cores. 
We expect that, to a first approximation, angular momen- 
tum is conserved minutely during the early stages of col- 
lapse, and this leads to the formation of a protostellar disc 
and/or to fragmentation into a multiple system. Thus, un- 
derstanding the formation and evolution of protostellar discs 
is a critical stage linking the initial conditions within prestel- 
lar cores to the final outcome of the star formation process. 

In this paper, we simulate the collapse of isolated, 
rigidly rotating cores, using the 3D OpenMP-parallelise d 
SPH code VINE |Wetzstein et al.ll200Sl : iNelson et al.ll2008h . 
All cores have the same mass and the same in itial density 
profile, given by a Bonnor-Ebert sphere (BES) (|EbertJI 19551 : 
lBonnorlll956t ): they are distinguished only by their angular 
momentum. We study the influence of core angular momen- 
tum on the properties of the protostellar discs and proto- 
stars that condense out of such discs. In a companion paper 
(Walch et al., in prep.), we explore and contrast the prop- 
erties of protostellar discs and protostars condensing out of 
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turbulent cores having the same net angular momenta as 
the rigidly rotating cores considered here. 

Previous numerical work on the collapse and fragmen- 
tation of cores has demonstrated that the critical ingredients 
are the net amount of angular momentum and its distribu- 
tion, the initial density profile, and the treatment of ther- 
modyna mics (for a summary s ee [Goodwin et al. 2Q0( $). For 
examp le J Mivama et alj (| 19841 ) and lBurkert fc Bodenheimerl 
( 199 3) examine the evolution of cores with a uniform ini- 
tial density and an isothermal equation of state. They 
find that the outcome depends on the parameters a = 



J\U, 



GRAV 
ROT 



and P = E/rot/I^gravI. where U, 
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tional energies of the initial core. Specifically, iMivama et alj 
(1981) only find fragmentation for a[3 < 0.12. 

However, observations of prestellar cores suggest that 
they are centrally co ndensed, and we l l represented by a 
Bonnor-Ebert profile (|Alves et al.ll200ll : lAndre et al.ll2004l '). 
which is the equilibrium configuration for an isothermal 
cloud. Centrally condensed clouds appear to be more sta- 
ble against fragmentat i on than cores with a uniform den- 
sity llLvnden-Belllll964l : lBosslll993l:lBurkert fc Bodenheimerl 



Il99tj ). For a review see (Bodenheimer et all (2000) 

To date, the only parameter studies of the collapse 
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of Bon nor-Ebert spheres have been performed using grid 
codes. iMatsumoto fc Hanawal l|2003T l investigate the col- 
lapse of 3 M© critical Bonnor-Ebert spheres with the den- 
sity enhanced by 10%. They invoke a barotropic equation of 
state, and investigate different total angular momenta, dif- 
ferent rotation laws, and bar mode perturbations of different 
strengths. Their simulations use a nested grid with mirror 
symmetry about the equatorial plane, and fine resolution 
in the centre (the location of the primary protostar), but 
quite coarse resolution in the outer parts (the disc region), 
particularly in the vertical direction. They conclude that 
fragmentation requires fi c t FF > 0.05, where f2 c and i FF are 
the initial angular speed and freefall time at the centre of 
the core; this is equivale nt to a l o wer li mit on (3. The study 
of disc fragmentation by Nelson (2006) suggests that these 
resu lts may be com p romis ed by poor vertical resolution. 

iBaneriee et~aH |2004) investigate the collapse of a 
2.1 M core and a 168 Mq core, using the FLASH AMR 
code. The cores are initially Bonnor-Ebert spheres with the 
density enhanced everywhere by 10%, and then a 10% m — 2 
bar perturbation superimposed; they are initially in rigid 
rotation. In order to capture the thermodynamics more re- 
alistically th ey use the tabulated molecular-line cooling pre- 
scription of iNeufeld et all (ll995T ). They perform only one 
simulation of a low mass core, which has f2 c i FF = 0.2 and 
produces a bar; they suggest that this bar might subse- 
que ntly fragme n t if th e simulation were followed further. 

iBoss et al.l (2000) have demonstrated that core frag- 
mentation is highly sensitive to the treatmen t of th ermo- 
dynamics. However, although iKrumholz et al.l ([2007) have 
performed 3D AMR simulations of 100 M and 200 Mq col- 
lapsing cores with radiative transfer, there has as yet been 
no comprehensive parameter study of low-mass collapsing 
cores using a realistic treatment of the gas thermodynam- 
ics. This paper is an attempt to perform such a study. 

In Sections 2 and 3 we describe the initial conditions, 
the numerical code and the constitutive physics. In Sections 
4 and 5 we present in detail the evolution of cores with, re- 
spectively, low and high Q Q , and in Section 6 we describe 
how the evolution changes as fl Q is increased between these 
limits. In the Appendix we develop an analytic model, and 
derive a condition for fragmentation based on the initial pa- 
rameters of a core. In Section 7 this criterion is applied to a 
sample of observed cores, and it is shown that they are all 
stable against disc fragmentation. In Section 8 we summarise 
our main conclusions. 



molecular hydrogen at 28 K, and therefore the isothermal 
sound speed is a Q = 0.34 km s" 1 . The central density is 
p c — 10~ 18 gem -3 , and so the mass and radius are 



M = 1.1 



Rr, — 



G (4nGp c 



37j£ B = 17.000AU. 



(1) 
(2) 



Here ip(£) is the Isothermal Function and its first 

derivative (see IChandrasekhar fc W ares 1949); i/) B = if) (£ B ) 
and ijj' = ip' (£ B ). The core must be contained by an exter- 
nal pressure 

P EXT = l.lp c a„e - *B = k B 5.5 x 10 5 cm" 3 K, (3) 

delivered by a hot rarefied medium. The freefall time at 
the core centre (hence the timescale on which the central 
primary protostar starts to condense out) is £ FF = 67kyr. 
The thermal energy of the core is U THERM = 3M G a 2 /2, and 
the self-gravitational potential energy is 



U, 



1 GMl GMl 
■- 0.74 ■ 
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Hence the ratio of thermal to gravitational energy is 

0.74 ; (4) 



a 
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we note that a > 1/2, because of the external pressure term 
in the Virial Theorem. 



2.2 The initial velocity field 

The only core property which we vary is the initial angular 
speed, f2 . For a core in rigid rotation about the unit vector 
k, the initial velocity field is then 



v(r) = fi k A r . 
The moment of inertia of the core is 
J G = 0.28 M Q R 2 Q . 
Hence the total specific angular momentum is given by 



j D = 1.8 x 10 21 cm 2 s _1 ( 



io- 



(5) 



(6) 



(7) 



2 INITIAL CONDITIONS 

2.1 The initial core density profile 

The observed column-density profiles of isolated prestel- 
lar cores are often well matched by approximately critical 
Bonnor-Ebert spheres, i.e. equilibrium isothermal gas clouds 
which are close t o the critical state for gravitational conden - 
sation (see e.g. lAlves et all 1200 ll ; iHennebelle et al]|2003h . 
We have therefore modeled our initial cores as marginally 
supercritical Bonnor-Ebert spheres, truncated at dimension- 
less radius £ B = 6.9 (the critical value is £ B = 6.45) and 
then with their density increased everywhere by 10%, to 
ensure immediate contraction. The gas is assumed to be 



and the initial ratio of rotational to gravitational energy by 



°- 036 (10^1)' 



(8) 



Our choice of representative values of fl Q is informed 
by two constraints, (i) The total specific angular mo- 
menta in low-mass cor es are typically of order 10 21 cm 2 s _1 
l|Goodman et al J 1 19931 ). (ii) The non-thermal velocity dis- 
persions in low- mass cores are usually subsoni c, and at most 
trans onic (e.g. iBarranco fc Goodman! 1 19981 ; lAndre et al.l 
2007). Therefore the rotational velocity at the edge of the 
core should not be much greater than a a . The values of Q Q 
that we consider are listed in Table [l] the equator of the 
fastest spinning core rotates at Mach 1.5 . 
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Run 


n Q 


*o 


^28 




M D 


R D 


F 


io-^s- 1 


kyr 


kyr 


M 


M Q 


AU 


1 


0.6 


94.5 


114.8 


0.53 


1.17 


250 


N 


2 


0.8 


94.5 


115.4 


0.31 


1.39 


344 


N 


3 


1.0 


96.5 


118.4 


0.25 


1.45 


433 


N 


1 


1.2 


96.5 


119.2 


0.23 


1.47 


567 


N 


5 


1.5 


98.5 


123.3 


0.17 


1.53 


739 


Y 


6 


1.6 


102.4 


126.7 


0.24 


1.46 


1091 


Y 


7 


1.8 


104.4 


130.4 


0.16 


1.54 


952 


Y 



Table 1. Column 1 gives the ID of the simulation, and Column 
2 gives the initial rotational frequency, f2 . Column 3 gives the 
time at which the first protostellar object forms, t Q , and Column 
4 gives the time at which 28% of the mass is either in the pro- 
tostellar disc or in protostars, t 2S . Column 5 gives the mass of 
the primary protostar at this juncture, M*, and Columns 6 and 
7 give the mass, M D , and radius, H D , of the protostellar disc at 
the same time. Column 8 indicates whether the disc fragments. 



3 NUMERICAL METHOD AND 
CONSTITUTIVE PHYSICS 

3.1 The VINE code 

The simulations have been performed with the Tree- 
SPH co de VINE. A fu ll desc riptio n of VINE c an b e 
found in IWetzstein et~al1 (|2008h and iNelson et all (|200ST ). 
VINE is parallelised with OpenMP directives. It invokes 
a leapfrog integrator, and individual particle time steps 
with a CFL tolerance parameter of 0.1. Gravitational ac- 
celerations a re estimated using a tree, with opening angle 
9 = 0.005 (|Springel et all l200ll ). The gravitational soft- 
ening length equals th e hydrodynamical smoothing length 
(|Bate fc BurkertHl997T l, which is adapted so that each parti- 
cle has A/" NEIB = 50 ± 20 neighbours. Hydrodynamical forces 
are treated with periodic boundary conditions, but gravita- 
tional forces are not. Art ificial viscosity is treat e d usin g the 
standard prescription of iGingold fc Monaghanl dl98&f) with 
a AV = 1 and /3 AV = 2, plus the Balsara switch ( jBalsaral 
1 19951). We have experimen ted with time-varying viscosity 
( Morris fc Monaghanl Il997l ) . but find that it makes no sig- 
nificant difference. 



3.2 Equation of state and molecular line cooling 

For the purpose of calculating the gross thermodynamics, 
we assume that the gas is pure molecular hydrogen, with 
ratio of specific heats 7 = 7/5 and molecular weight p = 
2; thus for simplicity we are neglecting the fact that the 
rotational degrees of freedom of H 2 are not fully excited at 
low temperatures, and we are neglecting the contribution of 
helium to p. The equation of state is 



P = ( 7 -l)«, 

and the isothermal sound speed is 



pm p 



1/2 



(9) 



(10) 



In solving the energy equation we adopt the following 



with the constraint that the temperature is not allowed 
to fall below T = 9K. At high densities, p > p CRIT , we 
switch off th e radiative cooling, and the gas then evolves 
adiabatically. iNeufeld et al. I i|l995T) have computed the equi- 
librium abundances of key molecules and atoms (CO, H 2 0, 
H 2 , HC1, 2 , C and O), and the resulting net cooling rate 
due to line emission, A, as a function of density and tem- 
perature, using the local velocity gradient method to treat 
optical-depth effects. Neufeld has kindly supplied these rates 
in the form of a look-up table for densities in the range 
10 3 «S (n(H 2 )/ cm- 3 ) < 10 
10 < (T/K) < 2500. Most 
used a barotropic equation of st ate and therefore does not 
treat thermal inertia effects (e.g. Goodwin et al. 2004; Bate 



10 and temperatures in the range 
previous work in this field has 



---S ( e -K-^ 

199Sl;|Matsumoto fc Hanawal2 003), although lBaneriee et al 



( 2004) use the same procedure as us. By invoking line cooling 
we hope to capture more realistically the thermal behaviour 
of low-density material as it falls onto the protostellar disc 
and then spirals into the central protostar. Ideally radiative 
cooling should also take account of continuum emission from 
dust, in the regime where the gas and dust are thermally 
coupled, for example by using the algorithm developed by 
IStamatellos et a l. (2007). Unfortunately the computational 
machinery to treat both the line cooling regime and the dust 
cooling regime does not yet exist. We are currently develop- 
ing an algorithm to combine these regimes. 



3.3 Resolution 

The core is modeled with 430,000 SPH particles, 

each having m ass m SPH = 1.4 x 1O~ J M0. Following 

iBate fc Burkertl (|l997l ). the minimum resolvable mass is 
therefore 2A/" NEIB m SPH = 1.4 x 1O -3 M . For compari- 
son the minimum possible Jeans mass (corresponding to 
P = Pcrit and T = 9K ) is !- 5 x 1O" 3 M . Therefore the 
Jeans mass is always resolved. In addition there are 24,000 
ambient particles exerting the external pressure that con- 
tains the core (see Eqn.[5J. These particles also have mass 
m SPH = 1.4 x 10 -5 Mq , but they have much higher tempera- 
ture, and their density is adjusted so that they fill the space 
between the core and the periodic boundaries. We do not use 
sink particles. Instead we limit the CPU time by imposing a 
minimum smoothing length, h Mm — 2 AU. Consequently we 
are able to resolve the details of the protostellar disc, both 
radially and vertically, but not the central protostar. 



3.4 Definitions 



particle whose density rises above 
3 is presumed to be part of a protostar. 



procedure. At low densities, p < p, 



1 ■a _Q 

gem , we 



use the cooling rates computed by INeufeld et all l|l995h . 



Any SPH 

10 _11 gcm~ 3 is presumed to be part of a protostar. The 
primary protostar is the one that forms first, at the centre 
of the core, and t Q is its time of formation. Values of t Q 
are given in Table [T] t Q increases monotonically with Q Q , 
because more rapidly rotating cores take longer to reach the 
critical density, p* . We define the disc to be all material with 
density in the range 10 -16 gem -3 ^ p ^ 10 -11 gcm~ 3 , since 
in general this material has S> |tv| (where v<f, and v r are 
the azimuthal and radial components of velocity). 

The evolution of a core can usefully be divided into two 
distinct phases. During the Isothermal Collapse Phase ra- 
diative cooling is very efficient and the gas temperature is 
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100 1000 10000 

Radius [AU] 

Figure 1. Radial density profiles during the Isothermal Phase 
of Run 1. Reading from the bottom, profiles are shown at t = 
0, 41, 51, 61, 71, 81 & 91kyr. The axes are logarithmic, and the 
dotted line gives a slope of — 2. 



roughly constant. Throughout most of the isothermal phase 
centrifugal acceleration is a secondary effect, and the core 
collapses on an approximately freefall time scale. The Proto- 
stellar Phase starts as soon as a disc forms. The pattern of 
evolution during the Protostellar Phase depends critically 
on Q Q . In cores with high £7 , (a) formation of the pri- 
mary protostar is delayed until a substantial disc has built 
up, and (b) the outer disc is prone to fragment producing 
secondary protstars. Since the Protostellar Phase is shorter 
than the Isothermal Collapse Phase, we give times in the 
form t — t + At, so that At is time measured from the 
moment when the primary protostar forms. 

Simulations are compared at time t 2S , which is defined 
to be when 28% of the core mass (i.e. ~ 1.7Mq) has density 
p > 1CP 16 gem -3 , i.e. it is either in the disc or in a proto- 
star. 28% is the typical efficiency of star formation observed 
in cores (e.g. iMotte et alJll998t lAlves et al]|2007h . In order 
to follow the simulations further, we would need to relax the 
imposition of adiabaticity at densities p > 10 -13 gcm -3 , so 
that the energy dissipation resulting from the redistribution 
of angular momentum can be radiated away. Without this 
radiation the long-term viscous evolution of the disc is in- 
hibited. 



4 A LOW ANGULAR MOMENTUM CORE 

In this section we describe the evolution of a low angular 
momentum core, having S1 Q = 0.6 x 10 -13 s -1 , and hence 
j Q = 1.1 x 10 21 cm 2 s -1 and p = 0.014. This is Run 1. 

The Isothermal Collapse Phase. Fig. 1 shows the 
density profile at various times during the Isothermal Col- 
lapse Phase of Run 1. The collapse is roughly spherically 
symmetric (at least until the late stages) and approximately 
self-similar: the central region with a flat density profile 
shrinks, both in radial extent and in mass. The infalling en- 
velope with p cx | r | 2 extends ever further inwards. Density 
profiles are computed by arranging the particles in order of 
increasing [r|, then averaging the density of contiguous sets 
of 100 particles; evidently this works only as long as depar- 
tures from spherical symmetry are small. 



-17.0 -16.0 -15.0 -14.0 -13.0 -12.0 




-100 100 -100 100 



Figure 2. The density on 400 AU X 400 AU slices through the 
centre of the core during Run 1. The left-hand (respectively right- 
hand) column shows slices at ,2 = (respectively y = 0). The times 
shown are (a) t Q — 1.9 kyr, (b) t a + 6.2 kyr, and (c) t Q + 18.2 kyr. 
The density scale ranges from 10" 18 to 10" 13 gem , and con- 
tours are plotted at 10 -16 , 3 X 10 -16 and 10 -15 gem -3 . Velocity 
vectors are superimposed, and the insets show a lkms" 1 vector. 

The Protostellar Phase. Fig. 2 shows velocity vec- 
tors superimposed on false-colour images of the density field 
on 400 AU x 400 AU slices through the centre of the core 
during Run 1. These images correspond to times from just 
before the formation of the primary protostar at t = t Q 
to t — t 28 . They show the formation of the primary proto- 
star and its attendant protostellar disc. The density field is 
computed on a 64 x 64 grid by first evaluating the column- 
density, E, through a slice having thickness A = 10 AU, and 
then setting p = E/A. The velocity field is computed in an 
analogous way, but on a coarser 20 x 20 grid. 

Density. Fig. 3(a) shows radial density profiles on the 
equatorial plane, p(r, z = 0), from Run 1. Most of the mass 
interior to 10 AU is in the primary protostar. Outside this 
the disc extends to beyond 100 AU, with p cx r~ 2,5 . After 
t = t a + 8.1 kyr, the edge of the disc is very sharp. 

Temperature. Fig. 3(b) shows radial temperature pro- 
files on the equatorial plane, T(r,z — 0), from Run 1. 
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Throughout most of the disc, the run of temperature on 
the midplane approximates to T oc r~ 0,6 . This is inter- 
mediate between the scaling expected for a thin, non-self- 
gravitating, steady-state Q-disc [which is heated by viscous 
dissipation, has a negative ve rtical temperature gradient , 
dT/d\z\ < 0, and T EFF oc r" 3/4 |Shakura fc Svunvaevl 19731) ] 
and the scaling inferred from observations of protoplane- 
tary discs around classical T Tauri stars [which are mainly 
heated by stellar irradiation, have a positive vertical tem- 
perature gradient, dT/ d\z\ > 0, and typically T EFF oc r -1 / 2 
jDullemond et al.ll200lh ]. Our simulations pertain to an ear- 
lier stage when the disc is self-gravitating and should be de- 
scribed as protostellar rather than protoplanetary. At this 
stage the infalling gas is heated impulsively when it encoun- 
ters the accretion shock at the boundary of the disc; evidence 
for shock heating can be seen in the temperature spikes at 
the outer edge of the disc, on Fig. 3(b). Since the rate of infall 
onto the disc exceeds the rate of accretion through the disc 
and onto the primary protostar, the heating at the accretion 
shock dominates the energy budget of the disc. Subsequently 
the gas is heated further by compression and viscous dissi- 
pation as it spirals inwards and on to the primary protostar, 
and therefore the vertical temperature gradient is negative, 
dT/d\z\ < 0. Consequently we anticipate that, as a disc ad- 
justs from being protostellar to being protoplanetary, the 
10 fim silicate feature will switch from absorption to emis- 
sion, as the sign of the vertical temperature gradient changes 
from negative to positive. 

Radial velocity. Fig. 3(c) shows radial velocity pro- 
files on the equatorial plane, normalised to the local isother- 
mal sound speed, v r (r, z = 0)/a(r, z — 0), from Run 1. Once 
the disc is well established (i.e. by t — t Q + 8.1 kyr), the 
gas in the equatorial plane is in freefall, until it hits the 
accretion shock at the edge of the disc (100 to 200 AU) at 
~ Mach4. The gas away from the equatorial plane hits the 
accretion shocks bounding the disc above and below even 
faster, because these shocks are even deeper in the gravita- 
tional potential well of the core. Thus the gas entering the 
disc is typically shock-heated to at least 100 K. The radial 
velocity of the gas in the disc is very tiny because the ef- 
fective viscosity arising from gravitational torques is small, 
and therefore the gas takes many orbits to migrate inwards. 

Azimuthal velocity. Fig. 3(d) shows azimuthal ve- 
locity profiles on the equatorial plane, v^r, z = 0), from Run 
1. Throughout the disc the gas is in approximate centrifu- 
gal balance; because the disc is massive and self-gravitating, 
this results in a profile with v<f>(r, z = 0) ocr _1//4 (significantly 
flatter than a purely Keplerian ti</>(r, z = 0) oc r~ 1//2 ). Out- 
side the disc, the infalling matter changes from its initial 
solid-body v^(r, z — 0) oc r to v$(r, z — 0) oc r _1 , as it falls 
inside ~ 6, 000 AU. 

Scale-height. At each radius r in the disc, we can 
define the vertical scale-height as 



z{r) 



E(r) 



p(r,a = 0) 



(11) 



where E(r) is the azimuthally averaged surface-density (as 
seen from z = ±oo) and p(r, z — 0) is the azimuthally av- 
eraged volume-density on the midplane (z = 0). Fig. 3(e) 
shows radial profiles of the disc aspect ratio, z(r)/r, from 
Run 1. The aspect ratio is almost constant, at a value 
z(r)/r ~ 0.55. Thus the disc is neither flat nor flaring; 
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Figure 4. The masses of the primary protostar, Af*(i) (magenta) 
and the protostellar disc, M D (t) (black), as a function of time, 
for Runs 1 and 7. 



rather the iso-density surfaces are approximately conical. 
Since the midplane density varies as p(r, 2 = 0) oc r _5//2 and 
the scale-height varies as z oc r, the surface-density of the 
disc satisfies E(r) oc r" 3/2 . [This i s very similar to what i s 
inferred for protoplanetary discs bv lDullemond et al.l (|200ll ). 
but our discs are much younger and more massive.] It fol- 
lows that the mass of the disc out to radius r is roughly 
M DISC (r) oc r 1//2 , and hence the equilibrium azimuthal ve- 



locity should approximate to 

1/2 

/ I r i / . . I j" I \ 



-1/4 



(12) 



as indeed it does (see Fig. 3(d)). 

Disc growth. Fig. 4 shows the growth of the pri- 
mary protostar, M*(t), and the growth of the disc, M D (t) 
The growth rates are roughly constant at M* ~ 2.2 x 
10" 5 M Q yr~ 1 and M D ~ 5.0 x 10~ 5 M© yr" 1 ; hence the 
disc grows faster than the p rimary protost ar. The disc is 
marginally Toomre unstable l|Toomrelll964l . see Fig. 7), to 
the extent that weak spiral arms develop and contribute 
to the inward transport of material through the disc and 
onto the primary protostar. However, the disc is not self- 
regulating, in the sense that the rate of transport through 
the disc does not adjust so that M* ~ M D . The disc is quite 
concentrated: typically more than 60% of the disc mass has 
density p > 10 _13 gcm -3 . Radiative cooling times in the 
disc are sufficiently long that the discs undergoes gravo- 
acoustic oscillations in the vertical direction, in response to 
the continuing infall of matter. 



5 A HIGH ANGULAR MOMENTUM CORE 

In this section we describe the evolution of a high angular 
momentum core, having Q Q — 1.8 x 10 _13 s _1 , and hence 
j Q = 3.3 x 10 21 cm 2 s" 1 and (5 = 0.13. This is Run 7. 

The Isothermal Collapse phase. The Isothermal 
Collapse Phase of Run 7 proceeds in the same way as for 
Run 1, except that, due to the higher angular momentum, 
departures from spherical symmetry occur sooner, and are 
more extensive. As a result, the Isothermal Collapse Phase 
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Figure 3. Profiles from Run 1 at times t +0.2 kyr (i.e. just after the 
primary protostar appears), t Q +4.2kyr, t Q +8.1kyr (i.e. the stage at 
which the disc becomes well established), t Q + 12.1 kyr, t Q + 16.0 kyr, 
and t Q + 19.9 kyr (the end of the simulation), (a) The radial density 
profile on the equatorial plane, p(r, z = 0). The axes are logarithmic, 
and a line with slope — 2.5 is shown to mimic the run of density in 
the disc between 10 and 200 AU. The key gives the line-styles used 
to represent the different times, (b) The radial temperature profile on 
the equatorial plane, T(r, z = 0). The axes are logarithmic, and a line 
with slope — 0.6 is shown to mimic the run of temperature in the disc 
between 10 and 200 AU. (c) The radial velocity profile on the equato- 
rial plane, v r (r, z = 0). The velocity axis is linear, and velocity is given 
in units of the local sound speed. The radius axis is logarithmic, (d) 
The azimuthal velocity profile on the equatorial plane, v^ir, z = 0). 
The axes are logarithmic, and lines with slopes of —1/2 and —1/4 
are included to show that in the disc <x r^ 1 / 4 , corresponding to 
non-Kcplcrian differential rotation. Note that this plot has a larger 
radial extent than the others, (e) The radial profile of the disc aspect 
ratio, z(r)/r. Both axes are linear. 



ends with the formation of a protostellar disc, and the pri- 
mary protostar does not form for another 15 kyr, by which 
stage the disc extends to 700 AU. 

The Protostellar Phase. Fig. 5 shows velocity vec- 
tors superimposed on false-colour images of the density field 
on 2000 AU x 2000 AU slices through the centre of the core 
during Run 7. These images correspond to times from just 
before the formation of the primary protostar at t = t Q to 
t — t 2S . They show the growth and fragmentation of the 
protostellar disc. Densities and velocities are computed as 



in Fig. 2. The features that distinguish this high-f2 run 
(Run 7) from the low-f2 one discussed in the preceding 
section (Run 1) are (a) that the primary protostar grows 
much more slowly and acquires most of its mass via the disc 
(rather than by direct infall), and (b) that the outer parts of 
the disc are able to fragment creating secondary protostars. 

Density. Fig. 6(a) shows azimuthally averaged radial 
density profiles on the equatorial plane, p(r, z = 0) from 
Run 7. Most of the mass interior to 10 AU is in the pri- 
mary protostar. Outside this the disc extends to ~ 100 AU 
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Figure 6. Profiles from Run 7 at times t Q (i.e. the moment at which 
the primary protostar forms), t Q + 5.9 kyr, t a +12.2 kyr, t Q +16.1 kyr, 
t Q +20.1kyr, t Q +24.0kyr, t Q +28.0kyr, and t Q +31.9kyr (the end 
of the simulation), (a) The radial density profile on the equatorial 
plane, p(r,z = 0). The axes are logarithmic, and lines with slopes 
— 1.5 and —2.0 arc shown to mimic the run of density in the disc 
from 10 to 100 AU, and from 100 to 1000 AU, respectively. The key 
gives the line-styles used to represent the different times, (b) The ra- 
dial temperature profile on the equatorial plane, T(r, z = 0). The axes 
are logarithmic, and a line with slope — 0.5 is shown to mimic the 
run of temperature in the disc between 10 and 1000 AU. (c) The ra- 
dial velocity profile on the equatorial plane, v r (r, z = 0). The velocity 
axis is linear, and velocity is given in units of the local sound speed. 
The radius axis is logarithmic, (d) The azimuthal velocity profile on 
the equatorial plane, v^fr, z = 0). The axes are logarithmic, and lines 
with slopes of —1/2 and —1/4 arc included to show that in the disc 
i>0 oc r -1 / 4 , corresponding to non-Keplerian differential rotation, (c) 
The radial profile of the disc aspect ratio, z(r)/r. The axes arc linear. 



with p oc r -1 ' 5 , and then from ~ 100 AU to ~ 1000 AU with 
p oc r~ 2 ' . The disc in Run 7 is far more extended than that 
in Run 1, and its density profile is significantly shallower. 
The spikes superimposed on the density profiles are due to 
spiral arms and fragments condensing out of the disc. 

Temperature. Fig. 6(b) shows radial temperature pro- 
files on the equatorial plane, T(r, z — 0) , from Run 7. 
Throughout most of the disc, the run of temperature on the 
midplane approximates to T oc r -0 ' 5 , which is very similar 
to Run 1. However, in Run 7 the temperatures are lower. 
There are two reasons for this. First, because the disc is 



more extended, the material hitting the accretion shock at 
the boundary of the disc has fallen less far into the gravi- 
tational potential well of the core and is therefore moving 
more slowly; it is therefore not heated so much in the accre- 
tion shock. Second, again because the disc is more extended 
and therefore has lower surface density, but also because it 
is cooler, its cooling radiation can escape more readily. This 
in turn is the main reason why the disc in Run 7 is more 
prone to fragmentation. Since the rate of infall onto the disc 
greatly exceeds the rate of accretion onto the primary pro- 
tostar (see below) , the heating at the accretion shock again 
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Figure 5. The density on slices through the centre of the core 
during Run 7. The left-hand (respectively right-hand) column 
shows slices at z = (respectively y = 0). The times shown are 
(a) t Q + 11.8 kyr; (b) t Q + 29.7 kyr; (c) t Q + 33.7 kyr; (d) t Q + 
33.7 kyr, zoomed to show the first long-lived fragment; (e) t Q + 
38.0 kyr; and (f) t Q +40.0 kyr. In all cases the slices are 2000 AU X 
2000 AU. The density scale ranges from 10 -17 to 10 _12 gcm~ 3 , 
and contours are plotted at 10~ 16 , 3 X 10~ 16 , and 10 -15 gem -3 . 
Velocity vectors are superimposed, and the insets show a 1 km s _1 



dominates the energy budget of the disc. Further heating by 
compression and viscous dissipation ensures that the verti- 
cal temperature gradient is negative, dT/d\z\ < 0. However, 
because the outer disc is very cool, we only expect to see 
the 10 fj,m feature in absorption from the inner disc. The 
spikes superimposed on the temperature profiles are due to 
fragments condensing out of the disc and being heated by 
compression. 

Radial velocity. Fig. 6(c) shows radial velocity pro- 
files on the equatorial plane, normalised to the local isother- 
mal sound speed, v r (r, z = 0)/a(r, z = 0), from Run 7. The 
gas in the equatorial plane is close to freefall until it hits the 
accretion shock at the outer edge of the disc. Because the 
disc is so extended, the material hitting the accretion shock 
is not moving very fast; for example, at the outer edge of 
the disc it is moving with v r < Mach2. Thus the gas en- 
tering the disc is only heated mildly. The radial velocity of 
the gas in the disc is small, although there are significant 
perturbations in the vicinities of fragments. 

Azimuthal velocity. Fig 6(d) shows azimuthal veloc- 
ity profiles on the equatorial plane, v^(r, z = 0), from Run 
7. Throughout the disc the gas is in approximate centrifugal 
balance, but because the gravitational field is dominated by 
the mass of the disc, rather than the mass of the primary 
protostar, the profile is significantly flatter than Keplerian, 
i.e. < (— dlnv^/dlnr) < 0.25, except for where it is per- 
turbed by a forming fragment. Outside the disc, the infalling 
matter changes from its initial solid-body v$(r, z = 0) ocr to 
v ( j 1 (r,z = 0)ocr~ , as it falls inside ~ 10,000AU. 

Scale-height. Fig. 6(e) shows radial profiles of the 
disc aspect ratio, z(r)/r, from Run 7. The aspect ratio is 
almost constant at a value z(r)/r ~ 0.6 within the inner 
disc (r < 200AU), and it falls off slightly at larger radii, 
where it is strongly perturbed by spiral arms and frag- 
ments. Again it is neither flat nor flaring; the iso-density 
surfaces are approximately conical. Since the midplane den- 
sity follows 1.5 < (— d\n p/dlnr) < 2.0 and the scale- 
height varies as z ex r, the surface-density of the disc fol- 
lows 0.5 < (— dlnE/dlnr) < 1.0, and the disc mass fol- 
lows 0.5 < (din M D1SC / dlnr) < 1.0. The azimuthal veloc- 
ity should approximate to ~ (GA/ DISC (r)/r) 1//2 , hence 
< (— din vj, /din r) < 0.25, as indeed it does (see Fig. 
6(d)), modulo the fluctuations due to forming fragments. 

Disc growth. Fig. 4 shows the growth of the primary 
protostar, M*(i), and the growth of the disc, M D (t), in Run 
7. Until the disc starts to fragment, and then a fragment 
merges with the primary protostar, the growth rates are 
roughly constant at M 4 ~ 0.7 x 10~ 5 M Q yr _1 and M D ~ 
4.7 x 10~ 5 M yr" 1 . Thus the disc not only forms before the 
primary protostar, but it also grows much faster. 

Fragmentation. Because the disc is massive, extended 
and cool, it becomes Toomre unstable (Q T < 1; see Fig. 7), 
and strong spiral arms develop. Where the cooling time in 
the dis c is shorter tha n the dynamical t imes, the disc frag- 
ments (jGammiel 1200 ll ; iRice et all l2003ft . The disc is only 
weakly concentrated: typically less than 20% of the disc mass 
has density p > 10 -13 gem -3 , and most of this dense gas is 
located in the spiral arms and fragments. 

Individual fragments. The first fragment condenses 
out at t ~ t Q + 30 kyr at ~ 200 AU. This fragment can be 
seen on Fig. 5(b) and Fig. 8. At its inception this fragment 
has mass ~ 0.05 Mq. Subsequently it grows in mass, but at 
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Figure 8. The density field on 140 AU X 140 AU slices through 
the centre of the first fragment to form during Run 7. The slice 
thickness is 6AU and the time shown is t Q + 33.7kyr. The den- 
sity scale ranges from 10" 17 to 10" 12 gem , and contours are 
plotted at 10~ 16 , 3 X 10~ 16 and 10~ 15 gcm~ 3 . Velocity vectors 
are superimposed, and the insets show a lkrns" 1 vector. 
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Figure 7. The Toomre parameter, Q T , calculated for the face- 
on projections of Run 1 at to + 18.3kyr (left plot), and Run 7 at 
to + 21.5kyr (right plot). 



the same time it spirals inwards and at t ~ t Q + 38 kyr it 
merges with the primary protostar. Two further fragments 
condense out at t ~ t Q + 40 kyr, at ~ 400 AU and ~ 600 AU 
(see fourth and fith panel of Fig. 5). The simulation is ter- 
minated soon after their formation, so it is unclear whether 
they also merge with the primary protostar, or whether they 
survive as secondary companions. 

Circum-fragmentary DISCS. Fig. 8 shows velocity 
vectors superimposed on false-colour images of the density 
field in 140 AU x 140 AU slices through the first fragment to 
form. Densities and velocities are computed as in Fig. 2, but 
using co-ordinates measured relative to the centre of mass 
of the fragment. We see that the fragment is surrounded by 
a well-defined, rotationally supported circum-fragmentary 
disc (CFD). Very similar discs are found around the other 
fragments. Typically they have radii between ~ 40 AU and 
~ 60 AU. 

Transport of angular momentum. Although most 
of the core mass infalls directly onto the disc, the spiral arms 



in the disc are effective at transporting angular momentum 
outwards, by means of gravitational torques, thereby allow- 
ing material to spiral inwards and accrete onto the primary 
protostar. Fragmentation further enhances the transport of 
angular momentum. This is because CFDs rotate in the 
same sense as the mother disc. Consequently, on the inside 
of a CFD (i.e. the side closest to the primary protostar) 
the shear between the material in the CFD and the mate- 
rial in the mother disc acts to slow down the material in 
the mother disc, thereby reducing its angular momentum. 
Conversely, on the outside of the CFD (i.e. further from the 
primary protostar) the shear between the material in the 
CFD and the material in the mother disc acts to speed up 
the material in the mother disc, thereby increasing its an- 
gular momentum. 

Merging. Efficient transport of angular momentum by 
the CFD may cause a fragment to spiral inwards and merge 
with the primary protostar. However, it is possible that the 
inspiral is being artificially accelerated by our imposition of 
adiabaticity at densities p > 10~ 13 gcm -3 , since this forces 
a fragment and its CFD to remain hot, extended and there- 
fore very dissipative, rather than allowing them to continue 
contracting towards protostellar densities. 

DISC mineralogy. Whereas in Run 7 the azimuthally 
averaged disc temperature never rises above ~ 200 K, in Run 
1 it rises above ~ 500 K, and locally it reaches ~ 800 K, suf- 
ficiently high to produce crystalline silicates. Thus we might 
expect crystalline silicates to be present in discs forming 
from low-fi cores, and hence to give rise to characteris- 
tic features in the SEDs of single protostars. Conversely, we 
might expect crystalline silicates to be rare in the cold, ex- 
tended discs forming from high-fi cores, and hence their 
characteristic features should be weak in - or even absent 
from - the SEDs of multiple protostars. The existence of 
a link between the amount of crystalline dust which can 
be found in protoplanetary discs and the initial conditions 
within the p arental molecular cloud core has already been 
suggested bv lDullemond et al.l (|2006l ). 
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Figure 9. Disc formation as a function of core angular momen- 
tum. For all seven runs, we show the density of the discs within 
slices of thickness 10AU about the disc midplance (z = 0) at tfi na i 
on the same scale. From top (Run 1) to bottom (Run 7) the core 
rotation Qq and thus the total initial core angular momentum is 
increasing. 



6 CORES WITH INTERMEDIATE ANGULAR 
MOMENTA 

Referring now to the complete ensemble of Runs (1 to 7), we 
note that, as Q Q increases, there is an essentially monotonic 
shift in the outcome; small departures from monotonicity at 
high values of Q Q can be attributed to the stochastic nature 
of fragmentation. Key parameters are presented in Table 1. 

For low Q Q , a substantial primary protostar is formed, 
surrounded by a compact, thick, hot disc; the disc is unable 
to fragment, and simply accretes onto the primary protostar. 
As Q Q is increased, the growth rate of the primary protostar 
decreases, not least because it acquires an increasing fraction 
of its mass via the disc (rather than by direct infall). At the 
same time the disc becomes more massive, more extended 
and cooler. Eventually, for D. Q > 1.35 x 10~ 13 s _1 (Runs 
5, 6 and 7), the disc becomes unstable against fragmenta- 
tion, i.e. it sati s fies both the To omre and Gammie conditions 
(|Toomre|[l96i ;[G ammi el 120011 ). We are unable to ascertain 
whether all the resulting fragments merge with the primary 
protostar, or whether some of them survive as secondary 
companions. 

As Q Q increases, the time, t Q , at which the primary 
protostar forms increases, and its subsequent growth rate, 
M t , decreases. Consequently, the overall pattern of accretion 
shifts from one in which the primary protostar forms first, 
and then the protostellar disc (low f2 ), to one in which the 
disc forms before the primary protostar (high f2 ). 

Fig. 7 shows the Toomre parameter, Q T , calculated 
for the face-on projections of Run 1 at to + 18.3kyr (left- 
hand plot), and Run 7 at to + 21.5kyr (right-hand plot). 
In the warm, compact disc from Run 1, the lowest val- 
ues of Q T are in the range 1 < Q T < 3, so that the 
disc develops weak spiral arms that transport angular mo- 
mentum by exerting gravitational torques, but it does not 
fragment (cf. Lin fc Prin glc 1990; Laughlin fe Bodenheimei] 



1 19941 ; iLodato fe Ricdl2004 iMeiia et al.ll2005h . In the cool, 
extended disc from Run 7, Q T falls below 1, leading to the 
growth of massive spiral arms and then fr agmentati o n. Thi s 
is consistent with recent findings o f Rafikov (2005); 
iMatzner &: Levin! (l2005lh I Whitworth & : Stamatellos 
(|200ej) ; IStamatellos et al.l (120071); iKratter et al. 
(|2008h : lRice &: Armitagej (120091 1 : IClarkel (|2009h . who 
suggest that in the inner disc (r < 50AU) fragmen- 
tation is highly unlikely, whereas it is possible at 
larger radii. 



7 COMPARISON WITH OBSERVATIONS 

In the Appendix we develop an analytic model for the early 
growth and fragmentation of a protostellar disc, in terms of 
its mass, M Q , its initial radius, R Q , and its initial angular 
rotation speed, Q Q . We show that cores with Q F < 1, where 



Qf = 



6.1M, 



17, 000AU 



(i 



35xl0- 13 ! 



(13) 



should be prone to disc fragmentation, whilst material is still 
infalling from the core envelope. We can apply this c r iterio n 
to the samples of cores obser ved b vlGoodman et all <ll993h . 
iBarranco fc Goodman! (|l998'), and ICaselli et all (|2002r i~for 
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Figure 10. The parameter Q F (see Eq n, I13[l is plotted agains t 
core mass, M n , for th e core sample s oflGood man et all l|l993t ). 
iBarranco fc Goodman] lll998ll . and ICaselli et al.l J2002ti . Cores 
with Q F < 1 are predicted to undergo disc fragmentation. 

which vales of M Q , R Q and il are listed. The resulting 
values of Q F are plotted against M Q on Fig. 10. 

Although a few cores have Q F ~ 1, none have Q F < 1. 
We infer that none of these cores is strongly predisposed 
to fragment whilst material is still infalling from the core 
envelope. The observed angular speeds are too small, and 
the strong heating at the accretion shock (where the infalling 
envelope material meets the disc) acts to stabilise the disc. 

This inference must be qualified as follows, (i) The core 
parameters used are very uncertain, in particular the core 
masses, (ii) Many of the cores appear already to have formed 
protostars (as evidenced by their association with IRAS 
sources), and even those that have not may already be in 
an advanced state of collapse. Therefore the observed values 
of R Q and Sl may not be "initial" in the sense we have de- 
fined them, (iii) Our simulations do not take account of the 
low levels of turbulence observed in such cores, (iv) We do 
not rule out the possibility that at a later stage, once infall 
ceases, some of these protostellar discs cool and fragment, 
particularly if more efficient cooling due to dust operates. 



8 CONCLUSIONS 

We have performed simulations of the collapse of rigidly ro- 
tating prestellar cores. Each core is initially a mildly su- 
percritical Bonnor-Ebert sphere (£ B = 6.9), with its den- 
sity increased by 10%, and contained by external pressure. 
The gas is pure molecular hydrogen with T = 28 K (hence 
isothermal sound speed 0.34 x 10 5 cms _1 ), the initial cen- 
tral density is p c = 10~ 18 gem -3 , and hence the core mass 
is M Q = 6.1 M , its initial radius is R Q = 17, 000 AU, and 
its external pressure is P EXT = k B 5.5 x 10 5 cm -3 K. The 
initial angular speed is varied from Q Q = 0.6 x 10~ 13 s -1 to 
n o = 1.8 x lO'^s" 1 (corresponding to 0.013 < f3 < 0.12). 

A slowly rotating core collapses to form a relatively 
massive primary protostar, and then a protostellar disc 
grows around this protostar. The disc is massive, but quite 
hot and compact, and so it does not fragment; it sim- 
ply accretes onto the primary protostar. There is some 



observational evidence for the existence of massive discs 



,( e -K- 


Eisner et alj|2005l; Rodriguez et al.ll2005l;lGreaves et all 


2008 


[iMann & WilliamsH2009f). However, the phase we have 



simulated is short lived and still deeply embedded in the 
parental core, so our discs would be rather hard to observe. 
We conjecture that such discs will show signatures of crys- 
talline silicates out to disc radii of several AU. 



A rapidly rotating core collapses initially to form a pro- 
tostellar disc, and then a protostar condenses out at the cen- 
tre. The protostar acquires almost all its mass via the disc, 
and therefore it grows quite slowly. The disc is extended, 
cool and thin, and tends to fragment. 

On the basis of an analytic model, we obtain a condition 
on the mass, initial radius and initial angular speed of a core 
(Eqn. [13]) such that its protostellar disc will fragment, whilst 
material is still infalling from the core envelope. Applying 
this condition to observed samples of cores we predict that 
such fragmentation is unlikely. 

In a companion paper (Walch et al., in preparation), 
we explore how these results are fundamentally changed if 
a prestellar core is turbulent, but has the same net angular 
momentum as the cores modeled here. 
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APPENDIX A: ANALYTIC MODEL 

We develop here a simple model to explain why rapidly ro- 
tating cores tend to fragment (and thereby presumably pro- 
duce multiple systems), whereas slowly rotating cores do 
not. For simplicity we assume that the initial core is cylin- 
drical, with height equal to its diameter and hence initial 
density p = M /2-tvRz. ; we make this assumption purely 
because it allows us to formulate analytically the structure 
of the resulting disc. The core is initially in solid-body rota- 
tion with angular speed £l Q , and we assume that during the 
Isothermal Collapse Phase, and during the early Protostellar 
Disc Phase, angular momentum is conserved minutely. Con- 
sequently the material initially in an annulus at radius r D 
has mass M ro — M o r 2 /R 2 "interior to it" in the disc, and 
ends up at radius r with angular speed f2 D (r) = 0, o r 2 o /r 2 . 
Centrifugal balance then requires that GM To /r 2 ~ Q. 2 D (r)r, 
and hence 



~ GM Q ■ 

The edge of the disc is therefore at 
3 4 r>2 



R 



R o n o 



D " GM Q ' 



and 



G 2 M 2 



It follows that the final surface density of the disc is 



En(r) _ M o 2-r dr 



G 2 M 3 



ttRI 2nrdr 2nR%Q.%r' 



(Al) 



(A2) 



(A3) 



(A4) 



Infalling material impinges on the disc at radius r with speed 



( GM To \ 
»d(0 - {—^— ) 



GM rn \ 2 GM n 



R%n Q 



We shall assume that it starts arriving at time 



GM r 



\GM 



(A5) 



(A6) 



and continues arriving at a constant rate until time 2t (r). 
Hence, in the time interval t D (r) < t < 2t D (r), the flux of 
material onto one side of the disc at radius r is 



E D (r) _ fGM c 



47Ti?2n4f5/4' 



2U0 V 
the instantaneous surface- density of the disc is 



(A7) 



E(r,t) = S D (r) 



t 



G 2 Ml J j_ _ _ 



2-KR%Q,%f \ f i 



t=—\ (A8) 



the ram pressure of the infalling material is 



%W^oW%(^ < GMoY M ° 



R 3 o J 4TvR o n 5 fi 



(A9) 
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and the flux of energy carried by the infalling material, and 
dissipated in the accretion shock, is 



It follows that at t = £„ 



0.29M o , and 



3>p(r>p(r) 



GM C 
"o 



(AlO) 



The disc cools mainly by molecular-line radiation. Under 
the conditions obtaining in the discs simulated here, the 
lines tend to be thermalised and optically thick. Hence the 
cooling rate per unit area of the disc, due to one species, M, 
is 

2 V 2 / M fc 9 3 T 9 (r,t) 



F M (r,t) 



&A u Z M h*<*H(r,t) 

2V 3 / M (fc B r)Xn 4 f 

3*A M Z M h*c e G 2 M 



t 

1 1 



(All) 



where I M is the moment of inertia of M, A M is the Einstein- 
A coefficient of the J = 1 — > transition of M, Z M is the 
fractional abundance of M (by mass). Since we are neglecting 
redistribution of angular momentum, we also neglect viscous 
heating of the disc; this is a reasonable assumption, given 
that the disc grows faster than the primary protostar (i.e. 
the viscous evolution of the disc proceeds on a longer time- 
scale than infall onto the disc). Hence the principal heating 
mechanism is dissipation of kinetic energy in the accretion 
shock, and therefore we estimate the temperature by setting 
F D (r,t) ~ F M (r,t). It follows that, during the period when 
the disc is growing, 



-FF ( Le - * = !)> M SYS 

since this is very close both to the juncture when the disc 
is most unstable, and to the juncture at which we compare 
simulations with different fi G , we will evaluate the stability 
of the discs at this time too. At t = t FF , the disc is complete 
out to J£ D /16, and is just starting to form at 7? D , so between 
these radii the disc is actively accreting. From Eqn. (|A17|) we 
see that the accreting disc is most unstable (Toomre Q T ~ 
1.07) at 7? D /16, and this is basically where the simulated 
discs with high S1 Q tend to fragment. The condition for the 
disc to be able to cool fast enough to fragment at R D /16 
(i.e. Gammie Q G < 1) then reduces to 

£l > 1.35 M<f* Ro^ ■ 



(A19) 

Eqn. ()A19|) agrees with the numerical results, which give 
fragmentation for Cl Q = 1.5 but not for Q Q = 1.2. It also 
indicates how the criterion for fragmentation scales with 
core mass and radius, and we exploit this to analyse ob- 
served cores in Section [7] The analytic model also accords 
well with the simulations in predicting that fragmentation 
should occur near J? D /16, which for those d iscs that frag- 
ment is between ~ 200 AU and ~ 700 AU (cf.lRafikovll2005l ; 
iMatzner fc Levinll2005l ; IWhitworth fc Stamatellosll2006l ). 



(k B T(r,t)f 



GM, 
Rl 



T 3 2 A M Z M h*c 6 M 2 



t 

T 1 



(A12) 



The Toomre and Gammie parameters for the actively ac- 
creting part of the disc are 

a(r,t)Q D (r) 



?tM) 



Q G (r,t) 



' 7rG£(r,i) ' 

AoolM) _ nGZ 2 (r,t)a(r,t) 



24F D (r) 



(A13) 
(A14) 



To proceed, we consider cooling by CO, for which 7 co 



1.4 x 10" 39 gcm 2 



10" 



and Z r . 



co — iu ° ' "CO 
we introduce dimensionless core parameters, 



6.1 M 



Ri 



R, 







17,000AU 



0.002: and 



10- 13 s- 1 



Then, 



T(r,t) ~ 60K i?~ 



M A - 



Q T {r,t) ~ 0.76 M ( 



19 . 29 . 4 

' 36 RZ 6 Mo r~ 



„ 1_ 

r 4 



(A15) 
(A16) 



Q G (r,t) ~ 0.19 mJ» R q ^ Cl~ 



23 



— -1 

r 1 



(A17) 



The mass of the system (primary protostar plus protostellar 
disc) is given by 



M SYS (t) 
M 




1 - J-f 4 

1 24 1 ' 



1; 

1 <i< 2; 

2 <t. 



(A18) 



